Phase diagram of force-induced DNA unzipping in exactly solvable models 
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The mechanical separation of the double helical DNA structure induced by forces pulling apart the 
two DNA strands ("unzipping") has been the subject of recent experiments. Analytical results are 
obtained within various models of interacting pairs of directed walks in the (1,1,..., 1) direction on 
the hypercubic lattice, and the phase diagram in the force-temperature plane is studied for a variety 
of cases. The scaling behaviour is determined at both the unzipping and the melting transition. We 
confirm the existence of a cold denaturation transition recently observed in numerical simulations: 
for a finite range of forces the system gets unzipped by decreasing the temperature. The existence 
of this transition is rigorously established for generic lattice and continuum space models. 

I. INTRODUCTION 

In recent years, micro and nanomanipulation of single biological macromolecules has become feasible due to a 
dramatic improvement of experimental techniques. By using devices such as optical tweezers [0,^), soft microneedles 
H and atomic force microscopes B, it is now possible to access physical and mechanical properties of fundamental 
biological objects, namely proteins, nucleic acids, and molecular motors, on the scale of individual molecules. Special 
effort has been devoted to the measurement of force-elongation characteristics of double stranded DNA molecules 
(dsDNA), determining its response to external forces and torques in the absence of enzymes. The mechanical unzipping 
of dsDNA structure by a force pulling the end of one of the two strands, the end of the other strand being anchored to 
some physical support, has been studied by Bockelmann et al. ElS], who have measured the average force along the 
opening of the two strands. Mechanical forces are in fact exerted on the DNA molecule by different enzymes during 
the process of DNA replication or transcription [^J^]. On the other hand, the double helical structure of dsDNA 
may be disrupted 'in vitro' by changing pH, solvent conditions and/or temperatu re p| . This transition is known 
as melting denaturation, and it has been long studied by theoretical physicists [|10|— [12|| . Instead, only very recently 
DNA mechanical denaturation has been the subject of theoretical studies. Most of these studies have considered 
so far a simple extension of the Poland-Sheraga model [jnj, in which the two DNA strands are homogeneous ideal 
polymer chains interacting with each other, introducing a constant force pulling apart the two strands by ac ting on 
one extremity of both strands. By using a mapping into a quantum mechanical problem, it has been shown (l^ [uj] 
that the opening of the two strands occurs only if the pulling force is increased to a critical value. The unzipping 
transition turns out to be a first order phase transition, whereas the melting transition is second order, in the ideal 
case, in d = 3 (recently both simulations and exact results have shown that the melting transition becomes first 
order, when considering mutually self-avoiding walks (SAWs) |l7],|l8|] ) . The heterogeneous case has also been studied 
with similar techniques p| . Very recently, Monte Carlo simulations of interacting pairs of self avoiding walks on 
a cubic lattice and of bead-and-string chains in the continuum space have determined the whole phase diagram in 
the force-temperature plane, revealing the existence of a re-entrant zipping-unzipping transition by decreasing the 
temperature for a finite range of forces . 

In this paper, we obtain exact analytical results for a class of simple lattice models of interacting pairs of homoge- 
neous directed self-avoiding walks. We extend the model introduced in ref. in D = 1 + 1 to the generic dimension 
case D = d + 1, and analyze both the scaling behavior of the first order unzipping transition and the multicritical 
scaling laws at the melting transition in the force-temperature plane. We also consider different versions of the model 
depending on whether the crossing of the two walks is allowed or not, or simply penalized by an entropy cost. In all 
cases, we find that the critical pulling force increases with temperature at low T, implying thus the existence of a cold 
unzipping transition |l9| . This seemingly paradoxical property is due to the competition between the energy gained 
by increasing the open portion of DNA and the entropy lost with the full stretching of the separated strands. Cold 
unzipping will be exactly proved for both ideal and self-avoiding chains in the lattice, and for a discrete chain with 
constant distance between consecutive beads in the continuum space. We also discuss the role of denatured bubbles 
forming in dsDNA opening, as opposed to the end-opening of the strands induced by the pulling force. By comparing 
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the phase diagram of the different models considered, we point out that the approximation of neglecting bubbles and 
considering only Y-shaped configurations shares many of the features of a mean-field type approximation, with an 
upper critical dimension d c = 4. 

Our paper is organized as follows. In Sec. II wc introduce the models of directed walks on the lattice in any 
dimension D and compute their thermodynamical properties, deriving the phase diagram in the force-temperature 
plane. In Sec. Ill we analyze the scaling laws at both the thermal melting and mechanical unzipping transition, and 
we highlight the features that are physically more interesting in our exact results (namely the existence of a re-entrant 
cold unzipping transition, our main result, and the role of forbidding the mutual crossing of the two strands). In 
Sec.IVa we discuss the physical meaning of re-entrance in our (lattice) models, generalizing the proof of its occurrence 
to the more realistic case of SAWs. In Sec.IVb we prove that also a discrete chain in the continuum space undergoes 
cold unzipping, and then argue why such an effect had not been observed in previous calculations |l4|-[l6| , performed 
in the limit of a continuum chain. Eventually, we discuss some related models (with generic SAWs) and give some 
additional technical details in the Appendices. 



II. EXACT SOLUTION OF THE LATTICE MODELS 
A. Introduction of the models and their behaviour for thermal melting 

We consider a simple class of models in a D — d + 1 dimensional hypercubic lattice Z D . The two strands of a 
homogeneous DNA molecule with N base pairs are mimicked by two SAWs, directed along the (1, . . . , 1) direction. 
The two chains have one end in common, while at the other end a force g is pulling in the (1,-1,0,..., 0) direction. 
The walks gain a binding energy — e (e > 0) every time the bases with the same monomer index (or same projection 
upon (1, . . . , 1)) interact, i.e. wrong base pairing is forbidden. 

The canonical partition function for two iV-step directed self avoiding walks (DSAWs) is: 

Zn{P,9) = ^2 PN(x,/3e)exp(Pg ■ x), (1) 
fez D 

where pN(x,(3e) is the canonical partition function of two directed interacting strands whose last base pairs are at 
relative distance x. The following recursion relation holds for ppf(x,(3e): 

D 

p N+ i(x,l3e) = ^2 Pn{x - ei + e ,l3e){l + (exp(/3e) - l)6 s j) , (2) 
*ii=i 

where e*j, i = 1, . . . , D, are the canonical euclidean versors of the D-dimensional space and 6 is the Kronecker delta. 

These and similar equations have been intensively studied, at g = 0, within simple models of DNA thermal melting 
|j~7].(l9]], and, in D = 2, within the context of random walk adsorption |20| and wetting problems pl| . Note that the 
choice of the (1, . . . , 1) direction along which the walks are directed, is crucial in allowing us to write local recursion 
relations. It is thus no surprise that the model belongs to the same universality class of random walks in d = D — 1 
dimensions p2[ . 

It is well known (sec |l7j for a recent review), that, in the absence of a pulling force, within the model defined 
by equation (0), dsDNA undergoes a phase transition between a low temperature double stranded phase and a 
high temperature denaturated state only for D > 3 (d > 2): in D = 2, 3 dsDNA remains double stranded at all 
temperatures. 

However, we note that in (j^) the two strands are allowed to cross each other, without any restriction. This seems 
rather unphysical, since every real chain has a finite "hard core" distance. Thus, one should consider also the case 
in which crossing between different strands is forbidden (recent studies of both homogeneous [ ^3| and heterogeneous 
p4j DNA melting transition have indeed considered the d = 1 case with forbidden crossing). Let us focus for the 
moment on the D = 2 case, with no pulling force, where the effect of forbidding crossing is most dramatic. Here, the 
direction (l/\/2, l/\/2) can be identified as "time", and its normal (l/\/2, — 1/\/2) as "space" (x). The model with 
crossing (w.c.) implies no restriction on the relative distance x, whereas the one without crossing (w.o.c.) implies that 
x cannot change sign, e.g. x > 0. The model w.o.c. is equivalent to surface adsorption models previously considered 
pofl . While the thermal melting of the two strands w.c. takes place at T c = oo, in the model w.o.c. T c = e/ log (4/3) 
as it can be deduced from ref. pCfl . To elucidate this point further, we have tackled an intermediate model, where 
we do not forbid crossing, but we make it disadvantageous, by assigning a cost V > each time the strands pass 
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through one another. With calculations similar to those reported below, we find that melting takes place at the 
critical temperature (as a function of V): 

-l 

(3) 

As expected, as V — > oo, equation (||) yields T c — > e/ log (4/3), whereas T c ~ e/V^ as V — > 0. 

To implement the non-crossing constraint in the model in generic dimension, we require that, if the DSAWs join, 
one coming from direction &i and the other from direction ej (with i ^ j), when they divide again, they cannot 
both proceed along their previous direction; i.e we forbid that the first walk goes along direction e*j and the second 
one along ej. This excludes one configuration out of D(D — 1) at the splitting point (the remaining D possibilities 
would lead to the zipping of the two strands). As regards the D = 3 thermal melting, again the model w.c. has 
T c = oo, while, if crossing is forbidden as described above, the critical temperature is T c = e/ log (9/8). This result 
was also found in a similar calculation by Rubin p5| ] . For D > 3 the models with and without crossing both undergo 
a denaturation transition at a finite, though different, critical temperature, so that the effect of forbidding crossing 
here is less important (the critical temperatures are the same at the leading order 1 / D in the D — > oo expansion, as 
expected). 
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B. Behaviour of the models at non zero force 



Let us now turn to the models with the pulling force g in generic dimension D, with g = (g, —g, 0, . . . , 0). We can 
find the asymptotic value of the canonical partition function ([j]) by locating jf7],[l9],^6) the singularity closest to the 
origin of its generating function: 

oo 

Z(z,0e,0g)=jT l z lf '£ i p N (g,0e)exp(p&.2), (4) 

iV=0 x 

which we will also refer to as the grand partition function. We stress that pn(x, (3e) is a generic partition function: 
its detailed form will depend on whether or not we allow crossing. 

To proceed, it is useful to partition the DNA molecule in ds helices (with the strands attached to each other) 
and bubbles, sequences in which the two chains share just the first and the last base pairs; we also single out the 
contribution of the unzipped end of the two DSAWs, i.e. the part from the last contact to the end. In this way, the 
grand partition sum (||) can be expressed as: 

Z(z,j3e,j3g)= — ^4^$ S(z,0g), (5) 

y ,h uy> I — Dzexp(3e) i cx p^ £ ) n(z\ w 

where we have defined: 

oo oo 

B(z) = J2 z N b N , S(z,f3g) = 1 + [ c *-i(x) ~ Civ-i(0)%] exp(/3g • x). (6) 

N=2 N=l x 

In (H), &at is the number of 27V-step bubbles, and cn(x) the partition function of two iV-step DSAWs never touching 
each other and having their last sites at a mutual distance x, and their initial sites at a relative distance ei —ej for some 
Iq 7^ jo- By summing over all possibilities for initial conditions (see equation (|^) below), b^ = ^Y^^j-i j=i c JV-2(ei — ej), 
so that we need to find an explicit expression for the cat's. The equations they obey are: 

D D 

c N+ i(x) = c N (x - ej + ej) - c N (6) ^ Ss^-g^ (7) 

Notice that in eq.(|^) acts as a sink or absorbing state, i.e. once the two walks join, they never leave. In this way, 
cn(x) with i/O counts the number of pairs of walks that never touch each other (and with the last monomers at a 
relative distance x), while cn(0) counts the number of remaining pairs of walks that at some point come into contact 
and then remain stuck together. This last quantity plays the role of an arbitrary constant and has to be subtracted 
away from the final counting as in (|^) (see p5| for another example in which this prescription was used in a D = 3 
example) . 
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In order to simplify the presentation from now on we restrict our calculations on models w.o.c. (detailed calulations 
for this case and the outline of calculations for the simpler models w.c. are deferred to Appendix A). The initial 
conditions for such models are: 

Co(x) = E 5 <H,3i-2j \ ~ S m -?j n i ( 8 ) 

where ei and ej are the directions forbidden as described in the previous subsection. 

By performing a Fourier and a discrete Laplace transform on eq.(J^), and by using the relation be twe en 6 at and cat 



given above, we can derive an explicit expression for the bubble generating function B(z) (see eq. (|A6| )). The singu- 
larities of the second denominator of the grand partition function (|5|) are necessary to find the melting temperature. 
We find that the equations locating these two singularities are : 

*i = JyH (9) 
exp(-/3e) -Dz 2 = B (z 2 ) . (10) 

The first singularity z\ is leading to the usual random walk behavior in the absence of any interaction. The second 
singularity z 2 (/3e) is a function of the strength of the attractive interaction between the two strands, determining thus 
the behaviour in the native zipped phase. At g = 0, the critical temperature for thermal melting is obtained when 
the two singularities above coincide {z\ ~ z 2 ). Unlike the models w.c, which have T c = oo in D < 3, those w.o.c. 
have a finite critical temperature in any D > 2 (see Appendix A for further details). 
As regards the force-dependent third factor in (^), it reads: 

S(z, 0g) = zS sing (z, 0g) + l- zc(z, 0). (11) 

A third new singularity, depending on the external force, and fundamental in our calculations arises when computing: 

d D q 



S S ing{z,(3g) = rp c (z, q) V] exppg - iq) ■ £\, (12) 

J[-k,-k] d y 27r > 7 



where c (z, q) and c(z, 0) are given in eqs. (A4,A5). Notice that in eq.(|l2|) one can immediately compute the integrals 



on dq3 . . . dqD, by using the well known identity 2t:6 (q) — ^2 X exp (— iqx), so that one is left with the double integral: 
77T\ l ■^xc(z,q 1 ,q 2 ,O,...,O)Y2exp[(0g-iqi)xi + (-0g-iq 2 )x 2 \. (13) 

-* (^) (2t) £% 2 

We now give details for the evaluation of the inner integral only, the outer one being equivalent. Defining h{q\, q 2 ) = 
(z, qi, q 2 , 0, . . . ,0)j , we aim at proving that the complex translation qi —>■ qi ~ i/3g can be performed in the 
integral /* w ^^zl^Ell, Not 

ice that h is a periodic function of q± . We extend qi to complex numbers and consider 
the contour 7 as shown in Fig. 1. By the residue theorem, one can write: 

n dq± exp (— iqiXx) i f^ 139 idy exp (— iirxi + X\y) 



2n h(q!,q 2 ) J Q 2tt h(-K + iy,q 2 
T dqi exp (— iq^Xi — Pgxi) | f idy exp (inxi + X\y) 
_ n 2ir h{qi - i(3g,q 2 ) 



f idy exp (iirxr + x x y) 
J -3a 2?r h{-n + iy,q2) ^ 



'09 — ' 9o 

where the sum runs over all poles qo of H = ^^fe^y 1 ^ ( & t fixed q 2 ) inside the contour 7 and Res(7i, go) is the 
residue of 7i at go- Note that due to the periodicity of h the second and fourth terms in the previous equation cancel 
and thus the complex translation can be performed provided that no pole of Ti. can be found inside the contour 7. 
An equivalent conclusion can be reached for the outer integral on dq 2 . The condition of having no poles inside the 
contours of integration is satisfied as long as z < z 3 ([3g), where: 

z 3 (Pg) = n (15) 

{D - 2f + 2 + 2 cosh (2pg) + 4 (D - 2) cosh {fig) 
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is found by solving h(—i(3g, if3g) oc (l—zf(—i(3g)) — 0, with f(q) defined in eq. (A.2). In this case, we find consistently 

that z 3 ((3g) is just the singularity in the resulting expression S(z, fig) = 1 + z c (z, —i0g) — c(z, 0) . 

As z 3 is always smaller than z\ (for g ^ 0), the singularity closest to the origin has to be determined between Z2, 
controlled by the attractive energy e, and z 3 , controlled by the pulling force g. If z 2 < z 3 , the DNA molecule is zipped, 
otherwise it is unzipped. The free energy per monomer / and the average distance between the two ends (projected 
onto g = (1, —1, 0, . . . , 0)), < x g >, are defined as: 

/= lim T lQg ^* PN (f ' ^ 6XP ^ ' ^ (16) 



JV^oo N 

df 

n^oo ag 



<x g >= lim <g-x>=—^ z N (17) 



These quantities read in the thermodynamic limit: 

/=Tlogz 2 (/? e ); ^^=0 g< 9 c{T,e) (18) 
/= T log z 3 [fig) ; - z 3 (fig) [4 sinh (20g) + 4 (D - 2) sinh (fig)] g > g c (T, e) (19) 



where the critical force g c (T, e) is found by imposing z 2 (f3e) = z 3 (f3g), as given in eqs. ( |10|Jl5| ). The above equations 
show the existence of a first order transition at g — g c (T, e), if g c (T, e) > 0. The phase diagram for the models w.c. 
can be found exactly in the same way (see Appendix A for some details on this). It is interesting to notice that the 
singularity z 3 (f3g) does not depend on whether or not we allow crossing, whereas z 2 in the w.c. case is different from 
the w.o.c. case. This is enough of course to make the whole phase diagram different for the two kinds of models, as 
will be discussed in section III. 

We stress here that the main interest in using directed walks is that they are a sub-class of SAWs in the same 
dimension. However, qualitatively similar results can be obtained also with simple random walks (RWs) |2q] , but in 
that case it would be not physically meaningful to forbid crossing as done above. 

Experimentally, however, it should not be hard to set up an experiment closely related to the calculation we have 
performed, by stretching the ds molecule in one given direction before applying the external force. 

III. SCALING LAWS FOR THERMAL MELTING AND UNZIPPING 

Let us focus for concreteness on the result for the model in D = 2 w.o.c, in order to analyze the scaling laws of the 
system. The phase diagram, as found also in [19], explicitly reads (see Fig. 2): 



g c (T,e) = | cosh^ 1 



1 



2 y/l - exp (-pe) - (1 - exp (-/3e)) 



(20) 



The model exhibits a first-order "unzipping" transition if we move at a fixed value of T < e/ log (4/3), as shown 
in jl8[). As is shown in Fig. 2, lim^^g g c (T, e) = | and g c (T,e) attains its maximum at T = Tm — 0.9e where 
9c(Tm) — 0.68e > |. The transition is second order at g = 0; as in jL7) we find that close to the melting point 

<n>= JL-logZN-r-ifirN^), (21) 

where < n > is the mean number of native contacts, with r = (T — T c ) /T c , and <p t = 1/2 for the crossover exponent 
at thermal melting in the D — 2-case. The scaling function f(x) behaves as usual 

f(x) —>Q asx — > +oo, 

f(x)~x for|x|«l, (22) 

f(x) ~ — Ixl 1 ^* forx — > — oo, 

such that < n >~ N^ 1 at the transition. 

The exact expressions for < n > and for < x g > can be found by inverse Laplace transforms of the quantities 
Z(z, fig, f3e), Z(z, f3g 7 (3e), g ® jg ^ Z(z, (3g, f3e). However, if we only need scaling relations in the thermodynamic 
limit, we can use the discrete Tauberian theorem (see p^] and Appendix B), which relates the critical behaviour of a 
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series to the asymptotic behaviour of its coefficients. By using this method, we find, in the vicinity of the unzipping 
mechanical transition: 



A i1 + A 2 t v ; 

9 A 1 -f + A 2 r v ; 

where Ai(T c , g c ) , A 2 (T C , g c ) are determined in such a way that d9c j^' e) = -^fft' T = _ 9 c) 1 9c, t = (T - T c ) /T c , 
and (f>i = 4>2 = 1 for the two crossover exponents, consistently with the fact that the unzipping transition is first-order 
(the explicit form of Ai(T c , g c ), A 2 (T C , g c ) is worked out in Appendix B). Note that Ai > 0, whereas A 2 is negative 
in the re-entrant part of the transition curve. The scaling functions hi t2 (x) behave in a similar way to f(x), and thus 
we get that < n >^ N^ 1 and < x g >~ N^ 2 are both extensive at the transition point. The physical interpretation 
is that a macroscopic portion of the chain is still in the double stranded state, but the rest of the chain is unzipped; 
these are just the two phases coexisting at the first order transition. This result will be used later to justify the use 
of Y-shaped configurations. 

It is instructive to compare this phase diagram with that of the same D — 2-case when we allow crossing. The 
transition line (plotted in Fig. 2), obtained from z 2 — z% (see Appendix A) is: 

g c (T, e) = Ttanir 1 [1 - exp (-/fe)] = ~ + ^ log [2 - exp (-/?e)] , (25) 

Both the models behave similarly near T = 0, namely they yield a transition to the cold denaturated state described 
in |lS[ ]. When we move at a constant force g such that g c (0) < g < Maxx [g c (T, e)], at low enough temperatures the 
molecule is unzipped, and it zips by increasing T. This feature of the phase diagrams seems at first sight paradoxical, 
since one would expect the critical force to decrease monotonically as the temperature is increased. The physical 
explanation of this result will be given below. In the model w.c, moreover, g c always increases, approaching e as 
T — > oo, and the two strands remain zipped for every temperature when g < e. In the model w.o.c, instead, the two 
strands unzip again by further increasing the temperature. As regards scaling, on the other hand, the two models are 
identical so that we can say that the effect of forbidding crossing is irrelevant in the renormalization group sense, but 
has a dramatic effect on the form of the critical line. 

Let us now discuss the results that we have obtained in higher dimensions (the critical lines for three and four 
dimensional models are shown in Fig. 3). Common to all models is the cold unzipping transition found in the 
D = 2 case. Moreover, the scaling laws at the unzipping transition do not depend on dimensionality in the crossover 
exponents. On the other hand, the thermal melting has dimensionality-dependent critical exponents (see pT| , p7| ). 
Consequently, the behaviour of the critical line g c (T, e) near the end point T = T c also depends on D. As T — > T~ 
the result that we find in generic dimension is: 

e/T-e/T c (D) D = 2,4 

,(T,e) I exp 



(e/T-e/T a (D)) 

[e/T - e/T c (D)f 2 [log (e/T - e/T^D))}- 1 / 2 D = 5 
[e/T - e/T c (D)] 1/2 D > 5 



D = 3 

(26) 



where a is a constant, and T c — oo in D = 2, 3 for models w.c. 



IV. PHYSICAL EXPLANATION OF RE-ENTRANCE 



A. Re-entrance in lattice models 



To obtain some physical insight into our analytical results, we compute the free energy of a Y-shaped molecule of 
DNA with a force pulling at the extremities (see Fig. 4). In this way, we neglect all the configurations with bubbles. 
Such an approximation is valid at low T and yields the exact result in the limit T — > 0. Since the configurations of the 
unzipped part are weighted by exp [fig ■ x), in this limit only the completely stretched configuration will contribute 
to the free energy for the unzipped part of the Y. In the limit of low temperature, the free energy is then: 



G 



F(m,N) (N ~ m)(e + T log fi)-2gm, (27) 

where m is the number of monomers in the unzipped part, N is the total number of bases and fi is the connective 
constant of a single walk. Note that eq. ( p7j ) (with g — g(l, 0, . . . , 0)) is valid also in the case of two self-avoiding 
walks (SAWs), since the connective constant of SAWs constrained to avoid the fully stretched unzipped part, is the 
same as for ordinary SAWs pi]j29] (for a more detailed discussion on SAWs we refer to Appendix C). By minimizing 



with respect to to, we find the critical force g c (T, e) such that (27) is minimum for to = when g < g c (T, e), and for 
m = N when g > g c {T, e): 

5c (T,e) T ^°i(e + Tlog M ). (28) 

This is indeed what we have found in our calculations, and confirms that re-entrance is a robust feature of lattice 
models, not depending on dimensionality or self-avoidance. In other words, ( p8| ) means that, at low T, it is more 
difficult to open a dsDNA helix as T is increased, because the energy gain obtained through the unzipping is more 
than compensated by the entropy loss, since there is only one possible completely stretched configuration versus fi N ~ m 
possibilities for the double stranded portion of the chain (we will see in the following section that the entropy loss in 
the continuum space exhibits a power law correction). For high enough T, on the other hand, also other configurations 
will contribute to the open portion, increasing its entropy, and the energy gain will eventually favour the unzipped 
state, except that in D = 2 w.c, where the presence of bubble enhances the entropy of the native portion at any T, 
as explained below. 

We have also calculated the phase diagram obtained by considering only the Y configuration: this amounts to 
putting the generating function for bubbles B(z) — in (^). As expected, the Y-approximation gives the exact 
behaviour rt2q) in the limit T — > 0, whereas it gets wronger and wronger as T is increased (see Fig. 2 where we plot 
the two-dimensional case explicitly given by g c (T, e) = ^ cosh -1 [exp (/3e) — 1] ). For the critical line near the melting 
point, this approximation yields g c (T,e) ~ (T c — T) 1 / 2 for all D, which from ( p6| ) is correct only for D > 5 (d > 4). 
This is because such an approach neglects bubbles, which are expected to be more and more relevant as d decreases 
(for d > d c — 4 once the two walks have splitted away, they basically never meet again) and T increases. In this view, 
it is clear that the puzzling behaviour found in the D = 2 model with crossing, where g c increases for all T, is due 
to the growing presence of bubbles: before DNA can be unzipped, we have to disentangle all the bubbles that are 
forming, which is harder and harder as T — > oo. 



B. Re-entrance in the continuum 

We now prove that also discrete chains in the continuum space undergo a cold denaturation for T — ► 0. Let us 
consider two iV-monomer chains in R d with no constraint on the directedness, with a force g pulling at the extremities 
and with constant unitary distance between consecutive monomers. As /3 — ► oo, bubbles can be neglected and only 
Y-shaped configurations contribute to the partition sum. The partition function for a Y-configuration then reads: 

N-l 

Z N (/3e, pg) = ]T Z* N _ m (/3e)ZZ m (Pg), (29) 

rn— 

where Z^_ m ((3e) = J Rd d d x\ . . . d d XN- m -i5 {\x\\ — l)...5 {\xN-m-i \ — 1) exp ((N — m)/3e) is the partition function 
of the zipped portion of the strands, and Z% m (f3g) = J Rd d d y x . . ■ d d y 2m S - 1) . ..6 (|y 2m | - 1) exp (/3g-J2i=i &) 
is that of the unzipped end. Chain discreteness is crucial in ensuring the validity of eq.(p9j) for f3 — > oo, since when 
one bubble has formed, the length of the stretched portion of the chain can be at most J2i=i~ Vi = 2 (to — 1) -A- . The 
integrals in (|2^) can be performed in any dimensions, and the evaluation of Z^ m (J5g) involves the modified Bessel 

function of the first kind of order 2^2. I n particular, Z^ m (fig) ~ ( f§ ) ex P (2 m A?)> an d there is a power-law 

correction with respect to the lattice result. We find that Zn ((3e, (3g) X)m=o ^~ m ~ 1 exp ((N — m)e)Z% m ([3g), 
where Qd is the surface area of the unit sphere in d dimensions. From this, the critical force is easily found: 



9c( r, e) ^°!-tW-!i„ g 



' : ;r- Tf ~ 



(30) 



where T is the Euler Gamma function. The critical line g c (T, e) increases at low T and re-entrance is present also in 
the continuum, even enhanced with respect to the lattice case. The leading term TlogT in equation (BG) (due to the 
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power-law correction in Z^ m {fig)) is indeed not present in the d — 1 case, when one correctly recovers equation ( psj ) 
for a lattice random walk with fi — 2. 

We finally wish to discuss the relation of the present treatment with previous work ( []l3|-[il| ) done on the unzipping 
of homo-DNA in the limit of a continuum chain. The ds molecule with the pulling force had been described by means 
of an effective hamiltonian, which, apart from an irrelevant center of mass term, reads: 

11= I ,ln\— (-) + V (r(n))-g.-\. CM) 




where f is the relative separation between the strands, b is the effective Kuhn length of single-stranded DNA and 
V is a realistic short range attractive potential, namely a delta function or a potential well (the two cases should be 
equivalent according to the standard quantum theory as long as J drV(r) is the same). The system described by (|3l| ) 
is equivalently represented by a quantum system with a non-hermitian hamiltonian. One finds [TiUl^ ] that there is a 
first order unzipping transition when the force reaches the critical value: 



4e (T)Td 

g c {T, e) = \l p , (32) 



where £o(T) is the ground state energy of the quantum hamiltonian obtained from ( pl|) when g = 0. In the quantum 
mechanical system, for g < g c (T,e) the ground state is a bound state, while for g > g c (T,e) the spectrum of ( |3l| ) is 
continuous. 

Let us focus on the d = 1 case, corresponding to our two-dimensional models without crossing constraints. It is 
well known that at g = both a symmetric square well and a delta function always have at least one bound state, 
meaning that DNA remains ds at all T. In Fig. 5 we show the critical force for the two potentials (obtained simply 
from (£52|)). It can be seen from the figure that at high T both the potentials saturate towards the same limit, as in 
our calculation for the model w.c. This is expected since the delta potential case can be recovered from our equation 

(||) in the limit of both continuum space and time (see Fig. 4). At low T, however, g c (T,e) ~ (s)^ 2 f° r a square 
well of width A, and is constant for the delta potential, so that re-entrance is present only in the former case, but 
with a behaviour different from that found in the lattice. This already signals that the low temperature behaviour of 
the solution obtained through the quantum mapping is rather unphysical. This is due to the fact that in the limit of 
continuum chain the chain constraint is modelled in a 'soft' way, by using a harmonic potential between consecutive 
beads along the chain. This interaction is usually assumed to be entropic (as in eq. (|3l|)) and thus effectively vanishes 
at T = 0. In that limit eq. ( pl| ) is describing a set of N 'free' particles moving in an effective potential determined by 
V and g. Consequently, the quantum mapping is expected to describe the system well except in the low temperature 
limit. Indeed, the re-entrance found for the square well potential is due to the unboundedness of the potential felt by 
the "free particles" as soon as g ^ 0. As T — > 0, fluctuations are not important, and the particles stay close to the 
minimum of such potential, which is — oo for a force whatever small but finite. This is enough to unzip the strands 

As Zhou has observed in [|16|, one can improve this model by placing a hard core either inside the potential well 
or at a finite distance from the delta (this is analogous in spirit to introducing the crossing constraint in our lattice 

models). Once again, one finds that for T near T c , which is now finite, both potentials predict g c (T, e) ~ (T c — T) 
as in the lattice model (see eq.(p(i|)), whereas for low T the delta shows no re-entrance and the potential well displays 
a T 1 / 2 behaviour. We note once again that the low temperature behavior of such 'quantum' models is an artifact, 
when considered as polymer chain models, due to the 'soft' enforcing of the chain constraint. In this respect, discrete 
chain models in the continuum space are more realistic, either with constant bond length or with harmonic springs 
between consecutive beads. In the first case we have indeed proved (see eq. (|30|)) the existence of cold mechanical 
denaturation. In the second case it has been shown to occur by means of numerical simulations nM . 



V. CONCLUSIONS 



To conclude, we have studied simple models for DNA mechanical unzipping induced by a pulling force. By using 
analytical techniques, the relevance of forbidding the crossing of the two strands and the role played by bubble 
formation in the denaturation process has been discussed throughout the whole force-temperature plane. The scaling 
properties of the system along the transition line have also been determined. The existence of a cold mechanical 
denaturation for a finite range of forces at low enough temperatures has been exactly proved in a general case for both 
lattice and continuum space models. The fundamental role of chain discreteness has been emphasized in comparison 
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with related models studied by means of quantum mapping techniques. Even if we neglected many important effects, 
such as base pairs heterogeneity, intrinsic helicity, wrong base pairing, and different stacking energy for the natured 
and denatured state, it would be interesting to experimentally test this prediction stemming from such a simple model. 

Furthermore, on the basis of preliminary results (analytical and numerical) we anticipate that even in the presence 
of quenched randomness in the contact potential our models display a re-entrant transition, so in this sense it is really 
a "universal" feature. 

Acknowledgements: We would like to thank F. Seno and S.M. Bhattacharjee for illuminating discussions. DM also 
acknowledges INFM funding. This work was supported by MURST grant. 
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APPENDIX A: COMPUTATION OF B(Z) IN MODELS OF DSAWS WITHOUT CROSSING AND OF 
THE PHASE DIAGRAM FOR MODELS WITH CROSSING 

We here first work out the details to obtain an explicit expression for B(z), defined in eq.Q. Starting from equations 
\j , we can perform a Fourier and a discrete Laplace transform to obtain: 







c ( z iQ) = Y zN ^l2 cxp(iq- x)c N (x) = ,_. (c^ (£) - zf (<f)c(z, 0)), (Al) 

JV=0 - Z A<7; 



where we have made the following definitions: 

D 



f (q) = 2 £ cos (<&-<?,) I f(q) = D + f (q) (A2) 

OC 

c(z,f) = z N c N (x); c Q (q) = ^exp(i<f- £)c (z) (A3) 

Using c(z, 0) = /r_ w ff ] D (27r)° c ( z ' ec l ua tion ( Al) is easily solved since, by permuting variables inside the resulting 
integrals, it is possible to see that they do not depend on the particular step which we forbid. We obtain: 

~, - 1 c (q) - zW 1 (z) [A(D)f (q)~ Cq (q)} 

cM = T^m TTzmW) ' (A4) 

c ^ = A ^iW^y (A5) 

where A(D) = 1 — is the reduction factor due to the crossing constraint. As a result, using B(z) = 



Y%,j-i,j=i c ( z ' ^ - ej), we obtain 



bM-aMtPwTTy (A6) 

1 + z 2 Vri(z) 



where 



The singularity Z\ — 1/ D 2 , leading to the usual random walk behavior, comes in when evaluating the above integrals; 
their denominator becomes negative as q — > for z > Z\. 

We give here also a brief outline of the calculation necessary to find out the phase diagram for the models of Section 
II of DSAWs when crossing is allowed. These models are simpler to solve than the corresponding models w.o.c, since 
it is not necessary to partition the molecule of DNA as done in the text in eq.(|5|), which allows us to write a more 
explicit formula for the grand partition function. Indeed we can simply start from the recursion relations (||) and, 
after Fourier and Laplace transform, find the expression: 

P(z,?,0e) = - l -— ) (A8) 

1 - zf(q) 1 - (1 - exp (-/3e)) W Q (z) 

for the quantity P (z, q, fie), defined as: 

00 

P (z, q, fie) = zN ex P *)PN& P e ) ( A9 ) 
n=o s 

with pn(x, fie) defined as in eq.(|l|) and 

f d D q 1 s 
Wo = / jr^. A10 

The final form of the grand partition function, for DSAWs w.c. in the presence of a pulling force, is: 

Z(z,fJe,(3g) = \ flUwn , (All) 

1 - z/z 3 ((3g) 1 - (1 - exp {-fit)) W {z) 



with z 3 {fig) defined as in eq.([L5j). Alternatively one can proceed exactly as in Section II, but with initial conditions 
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Co (a) 



D 

E 



(A12) 



The only resulting difference would be to have A(D) = 1 in eqs.( A4) and ( A6); after a bit of algebra one can convince 
hi mself that the partition function found by recollecting the different factors in eq.(|^) is the same as that found in 
eq.( All ). In doing this, one needs to compute the relations between the integrals Wq(z), Wi(z), W2(z), which may 
be done as, e.g., in |p5| . For example, the bubble generating function for both kind of models could also be simply 
expressed as 



B(z)=A(D)(l-Dz-l/W (z)). 



(A13) 



Changing the dimension dependent factor A(D) from A(D) = 1, when crossing is allowed, to A(D) = 1 — D ^_^ ■ 
when crossing is forbidden, is enough for the melting transition temperature to become finite for D < 3. 



APPENDIX B: TAUBERIAN THEOREM AND ITS APPLICATIONS 

The discrete Tauberian theorem (see Ej) states that the relations 



and 

cm N — -^r—L(N), (B2) 
(z*) T(p) 

are equivalent provided that p > 0, {cat} is a positive monotonic sequence and L is "slowly varying" in the sense 
specified in p7| . 

We apply this theorem here first to find the scaling laws of the unzipping transition to give an example of how the 
relations reported in section III can be recovered. When approaching the critical line at g = g c from the unzipped 
phase, the smallest singularity is z 3 (pg), as defined in eq. (Jig). The leading behavior as z — > z^ is: 

Z{z,Pe,pg) Z ^ 3 1 1 (B3) 

exp(-^e) - Dz 3 - B[z 3 ) 1 - z/z 3 



dZ(z,f3e,(3g) z^z~ 1 1 dz 3 /d(f3g) 

d{fig) ~ eM-M-Dz 3 -B(z 3 ) (l-z/z 3 f z 3 (f3g) ' 

dZ(z, /3e, (3g) z^z~ exp(— /3e) 1 

d(M ~ (exp(-/3e) - Dz 3 - B{z 3 )f 1 - zjz 3 ' 



(B4) 
(B5) 



where we have neglected inessential factors when g ~ g c . Consequently, the application of the Tauberian theorem 
yields: 



Z N ((3e,f3g) iV ~°° — — s n _ 57^*3"". ( B6 ) 



N—*oo 1 — jv 

exp(-/3e) - Dz 3 - B{z 3 ) ~ 

dZ N ((3e,0g) N ^ N _ N dz^[d{pg) 

d(Pg) ~ exp(-/3e)- Dz 3 -B(z 3 ) H z 3 (f3g) ' 
dZ N (f3e,f3g) n-^oo exp(-/3e) 



— — • ( B7 ) 
N . (B8) 



d(M ' (cM-fc)-Dz 3 -B(z 3 )) 2 3 

The canonical partition function Zn can be expressed in terms of previously defined quantities as Zn (/3e, f3g) = 
T.sPn{x, (3e) exp(/3<? ■ x) (see eqs.(|l|) and (|)). 

Note that at the critical line zi (/3e) = z 3 ((3g c ) , implying that the e-dependent denominator in the above equations 
vanishes, as the critical line is approached (that is 7 = (g — g c ) / g c <C 1 and r = (T — T c ) /T c < 1). In that limit the 
order parameters obey 
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A 2 (T C , g c ) = -/3 c g c [ D + — ^ c9 ^ \ z \ ^ c g c ) [ 4sin h (2/3 c g c ) +4{D-2) sinh {[3 c g c )} + (B12) 



dZ N /d((3e) n^oo exp(-/3e) 
Z N A X "f + A 2 t 

< Xg _ 9Z N /d(Pg) W^)^ 
Zjv z 3 (f3g c ) 

The coefficients Ai(T c , <? c ), A 2 (T C , g c ) determine the normal direction in the (log(T), log(g)) plane, and are reported 
below: 

Ai(T c , g c ) = (3 c g c (d + ^M^ll^ z | (/3c5c ) [ 4si nh (2/3 c5c ) + 4 (D - 2) sinh (/3 c5c )] , (Bll) 
/3 c eexp(-/3 c e) . 

As expected, it may be easily checked that, whatever direction we choose to approach the critical curve, the physical 
quantities as (B9) only depend on the projection of this direction on the normal to the critical line. To s how this, one 
only needs first to observe that criticality is achieved through the vanishing of the denominator in ( p39| ), and then to 
apply the implicit function theorem to this denominator. 

Recently, some authors ( p^ , [l4] |) have suggested, as an alternative order parameter for the unzipping transition, 
the number m of monomers from the last contact to the end. By using the same tools as before, we can find some 
quantities of interest such as the probability distribution of having m monomers "liberated" (in the equations below 
we suppose 1<< m <<N). In general: 

P{m) = £w«W,.,mPw (B13) 

where W are the pairs of directed walks, pw is their Boltzmann weight and n u . z . is the number of "liberated" monomers 
of the configuration. The possible enforcing of the crossing constraint does not affect the following equations, simply 
shifting the critical melting temperature T c . 
At zero force and T < T c , in the native state: 

P(m) oc exp [-mlog(zi/z 2 )]m a(d) , (B14) 

where m a ^ = m -1 / 2 in d =1, (log (m)) 1 in d =2 and constant for d >2. We recall that z\ = 1/D 2 , with D = d+1. 
At criticality: 

P(m) oc (TV — m)^- 1 m a M, (B15) 

where in the first factor there are logarithmic corrections in d = 2,4 to the power-law behavior controlled by the 
crossover thermal exponent <fit ||l7| . Above T c we have: 

P{m) oc (N - mf^m a{d) . (B16) 

with (3(d) = -|1 - d/2\ - 1 if d 7^2, in d =2 = l/(iVlog 2 (A)). At force different from zero, in the zipped phase, 

we get: 



while in the unzipped phase: 



if T < T r , and 



P(m) oc exp (-m log (z 3 /z 2 )), (B17) 
P(m) x exp (-(A - m) log (z 2 /z 3 )), (B18) 
P(m) x exp (—(N - m) log (zi/z 3 )), (B19) 



if T > T c . Just at criticality, the probability distribution is flat. 

For the Y-shaped configuration, we get, at force different from zero: 
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P(m) oc exp (—(N — m) log (z y /z 3 )) u.z. phase, (B20) 

P(m) oc exp (—to log [z^Jz y )) z. phase, (B21) 

P{m) oc const. at criticality, (B22) 

where z y = exp (— /3e) //x is the singularity controlling the free energy per monomer in the Y-zipped phase, and \i = D 
the connective constant of a single directed walk. At zero force we get: 

P[rn) oc exp (-(iV- to) log (z y /zi))m aW denatured phase, (B23) 

P(m) oc exp (— to log {z\ / z y ))m alyd ^ native phase, (B24) 

P(m) oc m a{d) at criticality. (B25) 



Let us give the details in a simple case, to get e.g. eqns . ( |B 2 0||3 2 l| , |B 2 2|) . One has: 

yPS 7/ree 

P( m ) = m (B26) 

Zn 

where is the partition function of two (N — m)-monomer chains with the end points in common (this is the 

original Poland-Sheraga model), Zf[ ee is the partition of two m-monomer chains with no base pair in contact, and 
Zn is given by eq. ([!]). We have, for the Y at non zero force, Z^ m — z y ^ , Z^ ee ~ z^ m , and 

Zn ~ z~ N native phase, (B27) 

Zn ~ z^ N above criticality. (B28) 



APPENDIX C: Y CONFIGURATIONS WITH SAWS 

In this Appendix we extend the analysis of section IV to a Y configuration where the two chains are generic SAWs 
and not directed walks. We take here g = g(l, 0, . . . , 0). 

Let us call c/v(0, x) the number of iV-step SAWs ending in x and starting in 0, and CN.g = J2s c Jv(0, x) exp (Jig ■ x) 
(so that cat = c N As a first step, we need to prove the existence of a connective constant /x^ for a single SAW 

in the presence of a pulling force i.e. cn,q ~ ^ for large N. The existence of a connective constant then follows J27j] 
from the inequality: 

cjy(0, x) exp (fig ■ x) < ^ c Nl (0, xi)c N ^ Nl (xi,x) exp (/3g- x{) exp {j3g ■ (x-xi)), (CI) 

which implies the subadditivity of cn From the connective constants of RWs and DSAWs, we can estabilish the 
following bounds (recall that here g = g(l, 0, . . . , 0)) for the connective constant of a SAW in d-dimensions: 

d-l + exp(/3|<f|) </Xg-<2(d-l) + 2cosh(/3|.g|). (C2) 

Let us now introduce the canonical partition function. It reads: 

JV-l 

Zn = 2J exp (/3(iV - TO)e)cjy- m -i ^ c 2m(^i.^) exp (fig ■ (x 2 - Xi)), (C3) 

TCI— Xi,X2 

where c' 2rn (xi,X2) counts the SAWs which do not cross the first N — m attached monomers of the Y and xi,X2 are 
the end-points of the strands. 

An upper bound is found if we allow the opened part of the Y to cross the first N — m zipped monomers: 

AT-l 

z up P er = ^ exp (/3(Y - m)e)cAr_ m _i ^ c 2m(xi,x 2 )exp(Pg- (x 2 -x\))- (C4) 



The grand canonical partition function defined from (C4) displays singularities in z = Z\ — -i? and in z = z 2 = 



KM 



where /i is the connective constant of a standard SAW in e?-dimensional space, i.e. fi = Mj = q. 
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To find a lower bound for Zjv, we restrict the Y configurations to those with the joined part of the Y constrained 
to lay in the half-space {x\x ■ ex > 0, e\ • A = 0} (the origin is in the bifurcation point), and the unzipped part forced 
to stay in the other half-space. 

N-l 

Zr=^exp(^-m)e)4ti £ ^ M (^,ft) exp ((3g- (x 2 - (C5) 

rn— X\ : X2 

where we have indicated with ( J£f t ( ri 9 ht ) number of n-step SAWs constrained to stay in the left(right) half 
space. As the connective constants can be rigorously proved to be the same even for walks constrained in such a way 
(see [ p7|p9[ | for details, it is important here that walks generated from a reflection through the hyperplane {x\x ■ e\ > 
0, e*i • A = 0} leave the scalar product g ■ x unchanged), the singularities of the grand partition function obtained 

from (pa) axe the same as those of X)jv=o Z^J' per z N . Thus the grand partition function singularities associated to 
Zn in eq. ( |C3[) are z = z\, Za given above. 

In summary, we have proved that a mechanical unzipping transition takes place when /i| = /iexp(/3e), similarly 
to the directed walk case. We suggest that a standard computation of (for instance through exact enumeration) 
could give an accurate phase diagram for the "Y approximation" with SAWs. 



The bounds in eq.(C2) give immediately two bounds within which the transition line g c (T,e) lies. In the T — > 
limit eq.(^8|) follows since both bounds give the same asymptotic transition line. This should be the same also for 
generic SAWs, since in the T — > limit it is sufficient to consider only Y-shaped configurations. This demonstrates 
the existence of cold unzipping at sufficiently low temperature for the generic self-avoiding case. 
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FIG. 1. Contour of integration 7 in the complex plane used for the evaluation of the integrals in (UA 
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FIG. 2. Phase diagram for two-dimensional models of DSAWs (e = 1 in the figure). The solid curve refers to the model 
with crossing, the long-dashed one to that without crossing and the dashed one to the model which considers only Y-shaped 
configurations. 
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FIG. 3. Phase diagrams in D=3 and 4 dimensions for models with and without crossing (e = 1 in the figure). It can be seen 
that the difference between the four dimensional models is negligible. In D=3 (model w.o.c), T c ~ 8.49e. Solid line: model 
D=3 w.c.; 

Long-dashed line: model D=3 w.o.c; 
Dashed line: model D=4 w.c; 
Dotted line: model D=4 w.o.c. 
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a) 



_ j 




b) 




FIG. 4. a)We show here an example of Y-shaped configuration for the DNA molecule: in this approximation bubbles are 
neglected. 

b)A completely stretched configuration for directed walks, which is dominant for T — > 0. 
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FIG. 5. Critical force in d — 1, found with the quantum mapping. The solid line is the result with a symmetric square well, 
the dashed one with a delta function. The parameters are chosen so that the integral of the potentials over space is the same. 
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